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We demonstrate that the probabihty distribution of the net number of electrons passing through 
a quantum system in a junction obeys a steady-state fluctuation theorem (FT) which can be tested 
experimentally by the full counting statistics (FCS) of electrons crossing the lead-system interface. 
The FCS is calculated using a many-body quantum master equation (QME) combined with a Li- 
ouville space generating function (GF) formalism. For a model of two coupled quantum dots, we 
show that the FT becomes valid for long binning times and provide an estimate for the finite-time 
deviations. We also demonstrate that the Mandel (or Fano) parameter associated with the incoming 
or outgoing electron transfers show subpoissonian (antibunching) statistics. 

PACS numbers: 



I. INTRODUCTION 

Various far-from-equilibrium relations, such as the Jarzynski relation [li [3, Q or the fluctuation theorem (FT) 
0, M, B 0) B B , obtained for classical systems during the past decade, provide new insights into the emergence 
of irreversible processes in physical systems [HI These relations follow from the observation that the ratio 

of the probability of a system forward and time-reversed trajectory is given by the exponential of a quantity, 
the trajectory entropy production, which when ensemble averaged, gives the entropy production in the system. 
These relations therefore quantify the probabilities of observing " non-thermodynamic" trajectories with decreased 
trajectory entropy production. These probabilities are infinitesimally small in the macroscopic world but are 
non-negligible in microscopic systems. However, the ensemble averaged dynamics always satisfies the second law 
(entropy production alwa ys g rows). With the recent progress in nano and mesoscopic sciences, these probabilities 
can now be measured [isl . Il4| . Because in the microscopic world quantum effects can be important, it is interesting 
to establish whether these fluctuation relations remain true in the quantum regime. This is still an open issue 
[isl [TgI . [TtI . [isL [igl . [20l . [2lj |. One of the major obstacles for a general formulation of a quantum FT is the lack of 
a clear concept of a measurable trajectory. It is therefore helpful to consider systems undergoing a well-defined 
measurement process. In the counting statistics of p hotons emitted by an atom or a molecule driven out of 
equilibrium by a laser field [12, [H, [H, [2^, [H, [S^, [H^ [H, [s^l, a trajectory picture is provided by the history of 
the detected photons. However, the reverse trajectory (where the laser mode absorbs a photon from the molecule) 
is not easily measurable. Electron counting statistics provides on the other hand a clear trajectory picture given 
by the history of the electron transfers between the system and the leads, where the reverse trajectory (electron 
moving against the bias) is a measurable quantity. Electron counting statistics in nanosystems has attracted recent 
interest fsil Isp . [H, [H, [H, [H, 113, [H, [H, 140| . l4l| . Individual electrons crossing quantum dots have been measured 
[43 , Id^ . |44|. |45| . Measuring the statistics of both forward and backward electron transfer events is essential to verify 
the FT and has recently been reported in Ref. [i^. Most studies have focused on the few lowest moments of the 
distribution. However, the FT is connected with the probabilities of large fluctuations which require the knowledge 
of the entire probability distribution. 

In this paper, we use the many-electron QME derived in Ref. I47II and the generating operator (GO) formalism in 
Liouville space developed for photon counting statistics (2^. l23l. [24 l25l [2^ [28j to calculate the FCS of electrons 



in biased quantum junctions. Our central formal result is an equation of motion for the GO whose solution can 
provide the full electron transfer probability distribution. Neglecting electron-electron interactions, this GO can be 
factorized into products of single orbital GO, each leading to a statistics similar to the one of the single resonant-level 
system studied in Ref. (4l|. By constructing the current GO from the full electron counting GO, we show that the 
probability distribution of the net number of electrons k entering the system from one of the system-lead interface 
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[k{t) = dTl(T), where /(r) is the current and t the measurement time also called binning time], satisfies: 

Pt{k)/Pt{—k) exp(/3eyA;), where eV is the difference between the left and right lead chemical potentials and 
(3 = {ki,T)^^ the inverse temperature. The probability of observing a current in the direction favored by the bias 
voltage V is exponentially larger than that of measuring the reverse current. For a non-biased junction the two 
currents are equiprobable. 

The QME presented in section [IT] is used to calculate the FCS of electrons in section IIIII In section IIVI we define 
the GO for the net number of electron transfer and show that the FT holds for long measurement times. In section 
|V]we derive closed expressions for the current and its power spectrum and for the Mandel parameter. In section IVTl 
we calculate the probability distribution for the net number of electron transfer for a model of two-coupled quantum 
dots and analyze the finite-time deviations to the FT. We also study the behavior of the average current and Mandel 
parameter as a function of the bias and temperature. Conclusions are drawn is section fVIII 



II. THE QUANTUM MASTER EQUATION 

The quantum junction is made of a system (e.g. quantum dot or single molecule) coupled to two leads. The system 
Hamiltonian reads Hg — ejcjcg, where cj (cj) is the Fermi creation (annihilation) operator for the s system orbital. 
The Hamiltonian of the left (right) lead is Hl = J2i ^I'^l'^i [Hr. = X^r ^rc\.Cr), where I (r) runs over all the left (right) 
lead orbitals. The entire junction Hamiltonian is H — Hg + Hl + -\- Ht where Ht — X^si/ \^suc\ci, + /i.e.] is the 
system-leads coupling [v — l,r). In [47|, we used second order perturbation theory in the lead-system interaction and 
projection operators in the number of electrons in the system to derive a QME describing the dynamics of the reduced 
system density matrix. The QME gives a Redfield equation |48i| in Fock space which only contains coherences between 
states with the same number of electrons. When the relaxation induced by the leads is much slower than the Bohr 
frequencies of the system, the fast oscillations in the interaction picture can be averaged out. This approximation, 
known as the rotating wave approximation (RWA) in quan tum o p tics, is often performed on the Redfield equation to 
guarantee that the final QME is of the Lindblad form [H, [H, Hfl [EH, [H . Our QME finally reads 

p" = -i[Hs, p''] + ^ {vssCsp"+^4 + WsscIp"~^Cs - yssclcsp"" -WssP'^Cscl+ h.c.) , (1) 

where p" is the reduced density matrix of the system projected into the n electron part of the Fock space. The 
complete system density matrix is given hy p = X^n^*"- When summed over n, Eq. ([T|) gives an equation for the 
total p which is of the Lindblad form [l^]. Vss^s and w^s's are related to lead correlation functions. Assuming a 
quasi-continuous spectra for the leads and neglecting the level shift contributions (which only modify the bare Bohr 
frequencies of the system) , we have 



where 

Jv) 



v^/ = ™,(e.)|ri^)(e.)p(l -/,(£,)) (3) 
wiyj = nny{e,)\T(yHes)\'fy{e.) • 

ny{e) is the density of state of the left or right lead (y = L, R) at energy e. fy{e) = l/(exp [f3{e — py)] + 1) denotes the 
Fermi distribution of the y lead and py is its chemical potentials. We assume pR = po and pL = po + sV , where V is 

the bias [see Fig[T]. w^^^ is the electron transfer rate from lead y to the s orbital and vi^^ is the rate for the reverse 
process. These obey the relation 



yiv) = e'3(^-M.)^(y) , (4) 



so that 



'^SS ^ss 

[L) (R) 



e'^^^ . (5) 
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III. GENERATING-FUNCTION FOR ELECTRON-COUNTING STATISTICS 

We consider a system with M orbitals and n spinless electrons, so that n = 0,1,..., M. The number of 
n-electron many-body states (hereafter denoted states) is given by C^^ = (^j^.jtil)\n\ ■ The total number of Fock space 

states is A'tot — '^ri=o^rf — 2*^- As a result of the weak lead-system coupling and infinite leads assumption, 
the Fock space coherences (FSC) between many-body states with different n are neglected and the number of 
elements of the full many- body density matrix reduces from Nl^ = 4*^ to N^^^ = 'EtLoiCnf- The space of 
the density matrices where FSC have been eliminated constitute our reduced Liouville space. By expanding the 
QME ^ in the eigenbasis of the system, the population dynamics obeys a birth and death master equation 
which is decoupled from the coherence dynamics. Electron-transfer events are counted by identifying the terms in 
the QME which are responsible for the transitions between the populations. Their sequence constitutes a "trajectory" . 

We shall recast the QME ([T]) in our reduced Liouville space as 

|p>- (i: + 7 + f)|p>= yW|p> . (6) 

Ai is the generator of the QME. C describes the isolated system dynamics 

C = -zJ^eMcs,-]- (7) 

s 

We denote the four possible processes depicted in Fig. [T] by 77 = 1, 2, 3, 4. 77 = 1 and 77 = 3 represent electron transfer 
from the system to the left and right lead whereas 77 = 2 and 77 = 4 represent the electron transfer from the left and 
right lead to the system. The orbital through which electron transfer occurs is denoted by s. F is responsible for 
electron transfers and is made of the non-diagonal terms of the generator which couple the populations 



where 1/ = (7;, s), = Ylt=i 'Etii^ 



f (1,,) = 27;(^)c, • 4 ; f(2^,) EE 2wi^^cl ■ Cs 

f(3,,) = 27;(f)c,-ct ; f(4,,) = 2ii;(f • c, . (9) 



7 describes the diagonal terms of the generator 



In analogy to F, we identify the various contributions to 7 from the different active orbitals 

l(l,s) = -~v[^s\c\csr]+ ; 7(2,5) = t«if'^[Cs4,-] + 

7(3,s) = -v[f[c\csr]+ ; l{i,s)=wif[cscl,-]+- (11) 

We can now calculate the full electron counting statistics using the formalism developed for photon counting 
statistics [H [H, S [H, [13, H H [131 . Starting with a trajectory picture of the QME evolution in terms of 
electron transfer histories, we shall calculate the probabilities of these trajectories and their associated generating 
function (GF). 

The system density matrix conditional to measuring k electron transfers during an interval of time t is denoted p'-''^ {t) . 
We use the compact notation defined in appendix [Al [see (jA7|) ]. k is a vector with components fc^. The probability to 
measure k electron transfers during a time interval t is obtained by tracing the conditional system density matrix 

P(i,k) =«/|/7('')(i)» . (12) 

The trace of an Hilbert space operator A, TrA, is denoted as a scalar product in Liouville space ^/|A:§>, where / is 
the unity operator. 
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The generating function (GF) associated with this probabihty distribution is defined as 

G{t,X) ^ 5]P(t,k)e^\ (13) 

k 

where A • k = \uku- Similarly, we define the generating operator (GO) as 

|G(t,A)» ^ ^|p('')(t)»e^'^. (14) 

k 

The GF is obtained by tracing the GO 

G(t,A) =</|G(t,A)> . (15) 
An evolution equation for the GO is derived in appendix El starting with the QME 

|G(<,A)>= W(A)|G(t,A)> , (16) 

where 

W(A) =M + Y^^e^" - ^)^- (17) 

is the generator of the GO evolution equation. Using the initial condition |G(0;A):^= |p(0)^, the solution of (jl6p . 
given in appendix [BJ provides the GF at all times 

G(i,A) =«/|e^W*|p(0)» . (18) 

The GF contains the entire information about the electron counting statistics. The probability distribution is obtained 
by inverting Eq. (|13p 

P{t,k)= dXG{t,iX)e-'^-^ . (19) 
Jo 

Moments of the distribution are given by derivatives of the GF 

gni Qn2 Qn 



= {KIK!---K:)t- (20) 

A=0 



We also define 



S{X) = - lim ilnG(t, A) . (21) 

t^oo t 



This will be useful to calculate the statistical properties of steady-state currents. 



IV. FLUCTUATION THEOREM FOR THE NET NUMBER OF ELECTRONS TRANSFERRED 

We will now focus on the statistical properties of the charge currents across the junction. We adopt the standard 
convention that the direction of the charge current is opposite to the electron transfers [see Fig.(IT])]. The number of 
electrons transferred via process rj through orbital s during a time interval t is given by 

1 

ku{t)^-- dTl,{T) , (22) 
e Jo 

where Iv{t) is the corresponding charge current \v — {rj, s)]. The net number of electron transfer events between the 
left (right) lead-system interface, through orbital s, during time t is 

*(L,s)(i) = fc(2,s)(i) - fc(i,s)(i) = -- / dTl^L^g-^ir) , fc(i^_s)(t) = A;(3_^)(i) - /c(4_s)(i) = -- / dTl^^R^^^^T) . (23) 

e Jo ^ Jo 
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where the charge current at the left [right] lead-system interface passing through the s'th system orbital is given by 
I(L,s){t) = I(2,s){t) — I{i,s){t) [I(R,s){'t) = I(3,s){t) ~ I{4:.s){t)]- The GF associated with the left and right net number 
of electron transfer G{t, Xl, Xr) is the GF of the FCS G{t,X) where Ai^s = —^2,s = Xl.s and X4^s — —X^^s = ^r,s- 
Defining the vectors (Xji) with components Al,s {Xr,s), we have 

G{t,XL,Xn) = ^ P(t,ki,kfl)e-^-''-e-^«"« . (24) 

For clarity, we hereafter consider the left GF. The right one may be calculated similarly. The GF for the left net 
number of electron transfer G{t, Xl) is defined from (f24|) by taking Xr — 0. Using Eq. (jB9p . we find that the generator 
yVs{X(L,s)) of the evolution of the single orbital GF Gs{t, X(l^s)) for the net number of electron transfer, via the single 
orbital s, through the left lead-system interface is given by 

Gs{t,XiL,s)) = c+(0)e^'+(^(-^=))*(g+i(A(i.,))+go+o(A(i,,))) 

+c_(0)e3-(^(--))*(gi^i(A(i^,))+go-o(A(i,,))) , (25) 

where, using (|B8[) and (O, the eigenvalues of the generator are given by 



(26) 



These possess the symmetry g{X(^L s)) — giP&V — A(i g)) so that, using Ss{X{l s)) = — limj^oo \ lnG's(i, X^^ s)) and 

Ss{X^L.s))^Ss{PeV -X(^L,s)) ■ (27) 

In appendix iBl we show that the many-body GF can be factorized into a product of single orbital GF [see (jBlOp ]. 
The GF for the total current (irrespective of the carrying orbitals) by setting A^^^-) = Xl in Xl) can therefore be 
written 

M 

G{t,XL) = ^e9'"(^-)* </|.g„,(AL)»<g™(Ai)|p(0)» , (28) 

m 

where gm{.XL), |3m(AL)^ and <C.gm(Ai)| are respectively the many-body eigenvalues, the right and the left eigenvector 
of the generator. Since the many-body eigenvalues corresponding to populations are made of 2^^ possible sums of 
single-body orbital eigenvalues ([26|l . they also satisfy the symmetry .g(Ai) = g{j3eV — Xl) so that, using (|2ip. 



S{Xl) = SiPeV - Xl) . (29) 

An important point is that the eigenvalues of the generator associated with the right GF are the same as those of the 
left GF so that S{Xl) = S{Xr). The eigenvectors will however be different and G(t, Xl) ^ G{t, Xr). This means that, 
in general, the electron transfer statistics at the left and right interface of the junction can be different. However, 
the statistical properties which can be obtained from S{X) are the same on the two interface. We will see that these 
include for example the FT [which follows from (|29|) ]. the steady state average current [see (|D2p ]. the current power 
spectrum at zero frequency [see (|D5|) ] or the asymptotic value of the Mandel parameter [see (jPSP ]. Hereafter, we 
therefore omit the L, R labeling in the corresponding quantities. 

In appendix [Cj we use the theory of large deviations to show that the symmetry (|29p implies at long times 

P{t,k) 

P{t, -k) ~ ^ ■ ^'^^> 

This is the FT for the net number of electrons k crossing the junction at each system-lead interface. Using ((27|) . we note 
that the FT also holds for the net number of charges which passed through each orbital k^y g) {y = L,R). A similar 
result was pointed out in Ref. [53| for a single resonant level in the large Coulomb repulsion limit excluding double 
occupancy. The FT implies that measuring electron transfers in the direction favored by the bias is exponentially more 
probable than the reverse process. The argument of the exponential is proportional to the nonequilibrium constrains 
of the junction PeV so that at equilibrium {V — 0) the two probabilities are identical. 
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V. AVERAGE CURRENT, MANDEL PARAMETER AND POWER SPECTRUM 

In appendix [D] we show how currents, moments and cumulants can be obtained from the GF. Using these results 
and the expressions for the eigenvalue with the smallest absolute value of the GF generators (|B8|) , we derive closed 
expressions for the steady state currents, their zero frequency power spectrum and the asymptotic value of the 
Mandel parameter. 



Using (jD2j) and (jBSp . the four steady state average currents through orbital s are given by 

(/(3.,)).t = -2e.(f , (/(,.,))3,^-2eu;(f)^^. (31) 

^ ' Vss +Wss ^ ' Vss + Wss 

The total current associated with the process 77 is obtained by summing over the orbitals (/^)st = X)f^(-^(i),s))st- The 
Fano parameter (t) and the closely related Mandel parameter M^{t) of each of these processes v = (77, s) are defined 
as 

F,(t)^NU(t) + l^^^^^^-^. (32) 

{ku)t 

The Mandel parameter vanishes for a Poisson process. Af < (M > 0) implies subpoissonian (superpoissonian) 
statistics. At steady state, using (ID7p . we find 

M,(oo) = iJ^^lls^ . (33) 

e Vss + Wss 

Another quantity of interest is the zero frequency power spectrum of the current Am, associated with the processes v 
[see Eq. (|D5p ]. Mi,{t) and Fi,{t) are easily related to it by A^^ = e{Iu)stFu{oo). The zero frequency power spectrum 
for the total current associated with the process rj is given by A^^ = ^w- The Mandel or Fano parameters cannot 
be expressed as such an orbital sum. It is therefore convenient to calculate them from S-q and {Iri)st- 
Similarly as for the t] processes, we can use (|26p to calculate the statistical properties of a given junction interface. 
Using (|D2[) . the average steady state current via the s orbital reads 

/ T \ It \ r\ I ^SS (^SS ^SS ^SS \ / j \ 

(/.).t = {IiR,s) )st = 2e -— . (34) 

\ ^ss ~r i^ss / 

Using (jD5p . the corresponding zero frequency power spectrum of the current is given by 

A _ \J-s/st 2(P- i j (35) 

Vss -I- Wss ' 

and A = Ef 



Wss V Vss + Ws 



VI. TWO COUPLED QUANTUM DOT MODEL 



We have calculated the probability distribution of the net number of electron transfer at the left lead-system 
interface for a model of two coupled quantum dots a and h. Quantities in local basis will be denoted by a tilde. The 
Hamiltonian of the dots in the local basis reads 

- ( t I ) (36) 
In the orbital eigenbasis, Hs is a diagonal matrix with eigenvalues 

= ^ ± sjC-^Y + n^ (37) 

The orbitals are labeled by s = 1,2. The two Hamiltonians are connected by a unitary transformation Hs = UHsU\ 
and similarly a = UaW and f3 — U(3W . We define the couplings between the leads and the dots in the local basis 
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and transform them to the orbital eigenbasis using U [see Fig. [T]. Since the two orbitals can be either empty or 
singly occupied, the system has 4 many-body states |00), |01), |10), |11). The many-body density matrix in the full 
Liouville space is thus a vector with 16 elements (4 populations and 12 coherences). In our reduced Liouville space 
it is a vector with 6 elements (4 populations and 2 coherences between |01) and |10)). The generator ([T7]) for this 
model is given in appendix lEl 

We have solved the GO equation for the total left current, trace the solution pS)) to get the GF, and finally 
calculate the probability distribution using (|19p . We assume that the measurement starts when the junction is at 
steady state. The parameters used in the numerical simulation are given in the legend of Fig. [Ij 

In Fig. [5^, we display the probability distribution of the net number of electron kj^ which have crossed the left 
lead-system interface for different measurement times and for a fixed temperature and bias. The logarithmic plot in 
Fig. [2)3 highlights the tails of the distribution. Positive (negative) kL represent electrons which move in the direction 
favored (unfavored) by the bias. As the measurement time increases, the average of the probability distribution 
moves to the positive direction linearly in time at a speed given by the steady state current. As expected, the longer 
the measurement time, the smaller the probability of observing a current flowing against the bias. Since the FT 
applies in the long times limit, this shows that the FT quantifies the rare fluctuations described by the tails of the 
probability distribution. 

Fig. shows the logarithm of P{t, kL)/P{t, —kL) for difi^erent measurement times. The longer the time, the closer 
the results from the FT. The numerical results suggest that for flnite times 



Pit,-k) 



g(/ieV-aOfc ^ (38) 



where limt^^oo ctt — 0. This form is not valid for very short measurement time, but seems to be a very good 
approximation for times after which the probability to measure at least a few electron transfers become signiflcant. 

To calculate at, we note that using ([24)l . Eq. p8|) implies that 

G{t, Xl) ~ G{t, (3eV - at - Xl) ■ (39) 
Since we do not consider very short times, (|28p can be very well approximated by 

\nGit,XL)^9moi>^L)t + B{XL) , (40) 
where toq is the index of the eigenvalue with the smallest absolute value and 

B{Xl) EE In ( «/|5mo(AL)»«5mo(AL)|p(0)» ) . (41) 

Since G(t, 0) = G{t, (3eV - at) = 1, we find that 

gm, {l3eV -at)t^ ~Bif3eV -at). (42) 

If we consider long enough times for which at <^ PeV and at/t « 0, using a first [zero] order expansion oi gmo{f3eV —at) 
[B{(3eV — at)] around [ieV and using (jD2p . we get 

^ B{PeV) 1 

at « . . - . (43) 

The average current can be calculated using (/)st — i-^s)st and ([M]) . Using B{Xl) — Bi^{L,s)): where 
^(A(L,s)) = In ( ^/|(ji_(A(L,s))3><C,g-(AL,s)|/5(0)» ) and p(0) correspond to steady state, we find that 



M 



B{f3eV) = ^ In . (44) 

[eP''^ vis wis + wis [vis + wis + wis )) 

Notice that since exp {i?(A(^ ^j)} > 1, then a* > (the equality only holds when /3eV = 0) 



Figure shows that our estimate for at is in excellent agreement with the values obtained by linearly fitting the 
results of Fig. [3^. It should be noted that the results presented on Fig. [312] correspond to small bias. For larger bias 
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the probability of the backward processes P{—k) becomes very small which limits the numerical accuracy. However, 
the GF is still numerically accessible for high bias. Fig. d] shows that the GF symmetry (|39p on which our methods 
relies is not perfectly preserved for larger bias. It is therefore expected that the accuracy of our method decreases 
with increasing bias. 

One can see on Fig. [2] that the probability distribution can be reasonably well fitted by a Gaussian. Deviations can 
be observed at very short times or for the tails of the distributions. The GF of a Gaussian probability distribution 
P{t,k) = (27rcr|)-i/2exp{-(fc - {k)ty/{2a^)} is given by G{t,\) = exp {A2ct|/2 - A(fc}t}. The nonzero solution of 
G{t,X) = 1 is Ao = 2(fc)f/<T^ 2/F{t), where F{t) is the time dependent Fano parameter associated with the net 
number of charge transferred through an interface. The GF has the symmetry G{t,X) — G(t, Aq — A). Using the 
results of section IVl to calculate F{oo), we find that for very long measurement times Aq = PeV + 0{V^). This 
indicates that it is only in the low bias limit that the Gaussian approximation can be trusted for describing the tails 
of the probability distribution which characterize the FT. This was confirmed by numerically studying the GF. We 
finally notice that our estimate of at is exact when the Gaussian approximation is satisfied, but can also be applied 
to non-Gaussian distributions as long as (f38|) or ([39|) remains a good assumption. 

The average steady-state current associated with the net number of electrons crossing a given junction interface is 
plotted as a function of the bias in Fig. [5^. This typical current- voltage characteristic shows that the current increases 
by steps each time the bias is large enough to make a new orbital contribute to the current (when eV = ei = 1.697 
or eV = €2 ~ 5.303). The current associated with the process i] — 1 and 77 = 2 is plotted on Fig. [5}d andO:. The 
current on Fig. [S^ is given by the difference between Fig. and Fig. [5}d. Temperature, when increased, has the 
effect of smoothing these steps and reducing the current because thermal fluctuations tend to equalize the forward 
and backward currents. Ohm's law is recovered for /3(e2 ~ £1) < 1- The zero frequency power spectrum associated to 
the total net current through a system interface is plotted as a function of the bias in Fig. [5ji. The Mandel parameter 
associated with the 77 = 1 and ry = 2 processes is plotted in Fig. ^ and [5]:. We see that the deviations from Poisson 
statistics are always strong and tend to vanish only when the currents associated with 77 = 1 and t] — 2 vanish. This 
indicates that the various type of electron transfer are highly correlated with each other. The Mandel parameter is 
always negative, indicating subpoissonian (antibunching) statistics. This has been experimentally observed in [46j . 

VII. DISCUSSION 

Many different types of fluctuation theorems (FT) have been derived for stochastic dynamical systems. These 
differ by the mechanism used to drive the system out of equilibrium. In the first case 0, H, [§], the system is 
closed and driven by a time-dependent force which makes the rate matrix of the birth and death master equation 
time dependent. When the driving stops, the system will eventually reach equilibrium because the rate matrix is 
detailed balanced 0] . In the second case 0, H, [§, HO] , the system is open and the the rate matrix is not detailed 
balanced. Even without driving, the system will eventually reach a nonequilibrium steady state. A third class of FT 
[5^ . [stI . [5^ considers the fluctuation of an entropy associated with the excess heat produced when a time depen- 
dent driving induces transitions between different nonequilibrium steady states. This paper focuses on the second case. 

So far, we have assumed that the two leads of the junction have the same temperature but different chemical 
potentials /i^ = /io + and fin = \xq. One could wonder what happens to the FT if one considers different 
temperatures for the two leads (3l = Po + = smd Pr = Po- In such case, the argument of the exponential on the 
r.h.s. of dn]) becomes x = PoAfi — Ap{es — Ho) + A/3A^. This imphes for the orbital GF that the analogue of the 
symmetry (P7)) becomes Ss{X(y^s)) = Ss{x — A(j,^s)). However, since x is different for each orbital (due to £5), the 
many body GF, which is given by the sum of the orbital GFs, does not possess the analogue of the symmetry (|29p . 
It is therefore only for a single orbital model that a FT P{t, k)/P{t, —k) = exp{a;/c} holds. 

In summary, we have applied the quantum master equation derived in (47| for calculating the counting statistics 
of electrons tunneling through a quantum junction made of a system embedded between two leads. Using a 
generating function formalism, we derived an evolution equation for the generating operator which allows to 
calculate the time dependent probability distribution of electron transfer events. This equation can be solved 
analytically because the many-body generating operator is a product of single orbital generating operators. We then 
demonstrated that the net number of electrons crossing a given system-lead interface satisfles a FT at long times. 
This implies that measuring a net number of electron transfer in the direction favored by the bias is exponentially 
more probable than measuring it in the opposite direction. Since the argument in the exponential is the work 
needed to transfer the measured electrons through the junction, this fluctuation theorem can be viewed as a Crooks 
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relation [3, y]. We furthermore described how the moments of the current distribution can be deduced from the 
electron counting statistics and gave analytical expressions for currents and power spectra. Numerically calcula- 
tions of the probability distribution for the electron counts for a model of two coupled quantum dots demonstrated 
that the FT becomes valid for long measurement times. A method to calculate the finite-time deviations was proposed. 

Several future extensions of this work are called for. The first one is to find if electron-electron interactions affect 
the FT. Another one is to investigate if the FT still holds in a system where the population dynamics couple to the 
coherence dynamics (e.g. a quantum junction externally driven by a laser). In analogy to the excess heat [svj . one 
could also consider fluctuations of excess currents, produced during transition between steady-states. 
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APPENDIX A: TRAJECTORY PICTURE FOR THE QME DYNAMICS 

Using the interaction representation, we can recast Eq. ([6]) as 

\pi{t)->=t{t)\pi{t):$> , (Ai) 

where 

|p,(i)»=Wo(0,i)|p(0» (A2) 

and 

f(t) =Wo(0,i)fWo(t,0) , (A3) 

where h(o{0,t) = exp [— A^o^]. A4o = JO + J describes the system dynamics in absence of electron transfer. The formal 
solution of (|Aip reads 

\pj{t)» = exp+ [ / dTt{T)]\pj{Q)»= H0» ' (A4) 

•^0 k=0 

where 

|pf(t)>- / dTk / 'dTfc_i... / 'dTif(rfe)f(rfc_i)...f(ri)|pKO)»= / dTf(r)|pf (r)> . (A5) 
Jo Jo Jo Jo 

Multiplying Eq. (|X5)) by ^^0(^,0), we get 

|pW(t)»= f dTk f 'dTk-1... / 'driWo(t,Tfc)mo(Tfc,Tfc_i)f...Wo(T2,ri)mo(Ti,0)|p(0)» . (A6) 
Jo Jo Jo 

This is the density matrix conditional to the transfer of k electrons between the system and the leads, irrespective of 
the type of transfer rj or orbital s. 

We shall denote the number of electron transfers of type 77 through the s orbital by where v = {r],s). The 
number of electron transfers of type 77 disregarding the orbital is fc^ = X)s=i ^{v-s) ^^'^ number of electron transfer 
through the orbital s disregarding the type of transfer is ks = X]»|=i ^(■n,s)- We have that k = k^. We define 

k = fc(l,2), • • • , fc(l,Af), fc(2,l), fc(2,2): • • • , fc(2,Af), ' ' ' , ^(4, 1) ' ^(4,2) , ' ' ' j ^(4,M) (A7) 

ks = fc(2,s), fc(3,s), ^(4,s) 

k,, = , fc(^,2) , • ■ • , k(ri,M) ■ 
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A trajectory {i',t)(^r) records (from left to right) the sequence of electron transfer events during a time r, by labeling 
them according to the type of transfer, the relevant orbital and the time at which a given transfer occurs (v, i){T) — 
(i^i, ti), (1^2, ^2), • • • ,{vk,tk)- A sequence {i^){t) is a trajectory where the transfer time is not recorded {v)(t) = 

We can now decompose pi{t) as 

k 

where the sum is over all the component fcj, of k and runs from zero to infnity. Using ([12]), the probability to measure 
k electrons transferred to the leads at time t is given by 

P{t,k) = ■^I\UQ{t,0)\p^^\t):^ . (A9) 
The master equation ^ preserves the trace so that the normalization condition P{t, k) = 1 is satisfied. 

We next define the elementary probability density of a given trajectory (v, T){t) which contains k electron transfer 
events at time ri , . . . , during a time interval t as 

n[(i.,T)(4)] =«Tr|Z:i'o(t,0)f,,(rfc)f,,_,(Tfc_i)...f,,(Ti)|pKO)» ■ (AlO) 
The probability of a sequence iv)(t) with k events is obtained by time integrating (jAlOp 

n[(i^)(t)] = / dTk f ' dTk-i ... r dnli[{v, r)(,)] . (All) 
Jo Jo Jo 

Using Eq. (|A6|) . and since our equation conserves the trace, we see that '}2{u)^t) ~ where 

sum over all possible electron transfer sequences (without keeping track of transfer times) . The trace of (|A5[) can now 

be written using l|A10[) and (jAlip as 

p(i,k)= n[(^)d' (A12) 

(!^)(t)6k 

where the summation is restricted to sequences such that the number of transfer events is k. 

Because a number k of electrons at time t can be realized by the four types of electron transfer processes (Fig. [Ij and 
via M different orbitals, we have 

\pf\t)^ = E f dTt,{T)\pf-^'^\T):^ , (A13) 

where li, = 0, • • • , 0, 1, 0, • • • ,0, where 1 is at the position v in the sequence. 
Using the interaction picture of the GO, we can rewrite (|A13[) as 

|G/(i,A)>= / dr(Ve^'^f,(T))|G/(T,A)> . (A14) 

By taking the time derivative and going back to the Schrodinger picture, we get Eq. ([TB|) . 

APPENDIX B: SOLUTION OF THE GF 

The generator of can be written as a sum of contributions of each orbital 

M 

W(A) = E^«(^«) ' (Bl) 

s=l 

where 

4 

W,(A,) =4+7s + E^^"'"'r(,,,) (B2) 
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and = -ie,[c|c,, •]. The GO can therefore be factorized as a tensor product of orbital GO 



M 



|G(i,A)>= J]|G,(t,A,)»® 



which evolve independently according to 

|G,(t,A,)»= W,(A,)|G,(t,A,)» . 
Projecting Eq. (jB4[) into the system eigenbasis and using the notation 

M 

<ni---nM;n'i---n^^|G(t,A)>= J]^ <n,; <|G4t, A^)>= ]J G„^;„^ (t, A^) 



(B3) 



(B4) 



(B5) 



s=l 



we get 



Go;o(i, As) 
Gi;o(i, As) 

VGo;l(i,As)/ 



>V(As 



Go;o(i, As) 
Gi;o(i, As) 
\Go;i(i,As) 



(B6) 



where 



>V(As 



\ 






"^^s ~ Vss ^ t^ss 





^ss "U^ss 



(B7) 



Since we work in the reduced Liouvillc space (where FSC are neglected), only the coherences such that 



J2^=i = Ssii "-S are kept in (fBS]) . 

The two eigenvalues of the generator (|B7I) corresponding to the population dynamics are given by 



± 



Vss + Wss 



/(A.) 



(B8) 



where 



/(As) = ^;(^){(e(^(l.'■)+^(^.-')) - l)^(^) + (e(^(i.=)+^(^.=)) - l)wif^ 

The two eigenvalues corresponding to coherences are obviously given by the two lower diagonal elements of (|B7p since 
population are decoupled from coherences. The right [left] eigenvectors of the generator can be easily evaluated. The 
two associated to the population read VVs(As)|(7±(As) g±{^s)\g±{^s)^ [^3±(As)|VVs(As) =^g±(As)|.g±(As)]. 
By tracing the solution of (jB6| we get the orbital GF 

(B9) 



Gs(t, As) -c+(0)e3+(^=)*(.9+i (As) +5o+o(A.)) +c-(0)ef-(^=)*(5i-i(As)+5o-;o(A.)) , 
where 5^^;„^(As) =<ns; ns|ff±(As)» and c±(0) =<5±(As)|Gs(0, As)>. The many-body GF ^ is given by 

M 

G(t,A) = []Gs(t,As). 



(BIO) 



This constitutes the solution of Eq. . If one does the spectral decomposition directly on the many-body generator 
one gets 



G{t,X) = ^ef-W* </|.g„(A)>«5„(A)|p(0)> , 



(Bll) 



where 5„i(A), \g„i{X)^ and ^g„i(A)| are respectively the many-body eigenvalues, right and left eigenvector of the 
generator. Each of these many-body eigenvalue is made from one of the 2^ possible ways of summing the orbital 
eigenvalues (IB8|) and the many-body left and right eigenvectors are tensor products of the single orbital eigenvectors. 
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APPENDIX C: FLUCTUATION THEOREM DERIVED FROM THE GENERATING FUNCTION 

SYMMETRY 

The following reasoning is based on the steady state FT for the entropy first obtained [gI, ^ and later extended for 
currents [sj . [sSj . The GF is associated with the probability distribution by 



Aft 



(CI) 



where we have introduced P(t, the probability that £^ = k/t takes a value in the interval [^,^ + d^]. The large 
deviation function (LDF) is defined as 



This definition follows from the ansatz 



where 



t^oo t 



(C2) 



(C3) 



lim -lnC(C,t) = 



(C4) 



We can then rewrite (ICID as 



G{t, A) = y d^Cit i)e-(«(«)+A?)t (C5) 

At long times, the main contribution to this integral comes from the value of ^, that maximizes the argument 
of the exponential. ^* is i 
integration, (|C5|) becomes 



of the exponential. ^* is therefore the value of ^ such that A = — -^1^=^.. At long times, using steepest descent 



Substituting (IC6| in ^ gives 



dCCi^,t)e 



2tt 



S{X) ^ R{C) + XC 



(C6) 



(C7) 



This shows that S{X) is the Legendre transform of the LDF. The LDF is given by the inverse Legendre transform of 
5(A) 



i?(0 = 5(A*)-A*e 



(C8) 



where ^ = ^|a=A'- Since 5(A) is convex downward [i.e. d^S/dX^ — — hmt_^oo((fc^)t — (^}t) < 0]i its Legendre 
transform is convex upwards. Using the symmetry Eq. (jCSp imphes that R{~C) — S{[3eV — A) + {[3eV — A)^, 
which together with i?(0 = 5(A) — Af leads to 



Substituting this in Eq. (|C3|) . we get 



ln-^./^ey^. + ln.^(^'*) 



Using den, this gives the FT dSO]) in the long time limit 



c(-C,t) 



(C9) 



(CIO) 
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APPENDIX D: CURRENT, POWER SPECTRUM AND THE GENERATING FUNCTION 

We show how to calculate the average currents, their zero frequency power spectra and their Mandel parameter 
from the GF. These results are used in section fVl 

To simplify the notation, we assume we have a probability distribution P(t, k) and its associated generating 
function G{t,\) — '^^P{t,'k)e'^'^, where the component of k and A are given by k^j and A^,. We also have 
S{X) — — limt^oo 7 lnG(i, A). The averaged number of charge rj is given by 



- f dT{I^{T))= ^G{t,\) 

e Jo dXn 



(Dl) 



A=0 



and the steady state current by 



St ^ lim - 

t^oo t 



dr{I^{T)) = e—S{\) 



x=o 



We also find that 



(D2) 



(D3) 



Since at steady state (/r,(Ti)/^' (r2))sf = (/,,(ti - T2)/^' (0))st, (/,,(Ti)}st = and (/^'(T2))st = (/,,')st, using the 

fact that 



_d d_ 



lnG{t, A) 



{kjjkjji)t {kri)t{kri')t 



A=0 



we find 



dXn dXr,' 



A=0 



dT(^(X„(T)V(0)),t-(/^),t(V),t 

dr{[I^{r) - {l.,)st][l,'{0) - {I,r)st])st 



(D4) 



(D5) 



Since the Fourier transform of the current correlation function is given by 



dre— ([/„(r) - [^(0) - t])st 



(D6) 



we see that A^^/ in (jPSp is the zero frequency power spectrum of the current correlation function A^j^y = A^j^y {lu ~ 0). 
This quantity is used to study shot noise [32] . 

The analogue of the Mandel parameter in photon counting statistic [l^, for the process y is given by 



The asymptotic value is given by 



A=0 



(krj)t 



dX 



-G{t,X) 



1 . 



(D7) 



A=0 



M^(oo) 





A=0 




A=0 



- 1 . 



(D8) 



For a Poisson process Ai^ = 0. The zero frequency power spectrum is related to the long time limit of the Mandel 
parameter by 



A 



e(/„),t(l + M„(oo)) 



(D9) 
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APPENDIX E: GENERATOR FOR THE TWO QUANTUM DOT MODEL 



We present the basic quantities needed in the reduced Liouville space to study the model of two quantum dot model 
presented in section HVl 

The orbital eigenbasis is denoted by {|ni,n2)}, where rii {n2) is the occupation number of the orbital s = 1 {s — 2). 
Using the notation (rii, n2|p|n'j, nj) = Pnin2;n[n'^, the density matrix in the reduced Liouville space is given by the 
vector p = (poo;00 , Poi:Oi , PiO:io , Pii;ii , PiO:Oi , Poi;io)^- The generator of our QME H]) for this model in this 
basis reads 



M 



V 



2(wii + W22) 


2V22 


2vn 








2W22 


-~2{V22 + Wii) 





2vn 





2wn 





-2{vii + W22) 


2V22 








2wii 


2W22 


-2{vii + V22) 

















-X 














-X* 



\ 



J 



(El) 



where X — vn + wn + V22 + W22 + — £2)- As expected, the populations are decoupled from the coherences. The 
generator is diagonal for the coherences and obeys a birth and death master equation in the population space. The 
generator for the GO evolution equation is given by W(A) = 



W22j 



( -2(u;n 


where the coherence part has been discarded since it is the same as for the generator of the QME 







-2(u22+wii) 2(e^(i.i)wi^^ +e-^(3.i)wif^) 

-2(wii+W22) 2(e^(i.2)42^ +e^(3.2)i;^2^) 

2(e^(2.i)w^^^ + e^f-i.DwJf ) 2(e^(2.2)u;^2' + e^'-^'^^wf^) -2(wii + W22) 
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FIG. 1: Schematic representation of the model of two quantum dots a and 6 coupled in series between two leads used for 

our numerical results. The upper part depicts the system in the local basis where the Hamiltonian of the dots reads Hs = 
"^^^Hijclcj, where Haa = Ea = 2, Hi,t = Et — 5 and Hat ~ H^^. = = 1- The coupling clement with the leads are 



Tal = 0.5 Tbr = 0.3 and Tar = Tu = 0. The lower part depicts the system in the cigcnbasis where the Hamiltonian becomes 
Hs = "^gtsclcs, with ei = 1.697 and £2 = 5.303. The coupling element with the leads become Tir = —0.479, T^i = 0.145, 
Tlr = 0.087, T2r = 0.287. We choose //o = and ny{es) = 7r~^ as well as e = 1, ?i = 1 and kb = 1, so that the units of energy 
and temperature is f2 and the time unit Q"^. The four possible types of electron transfer (and their associated currents) are 
represented by the big dashed arrow. 
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FIG. 2: a) Probability distribution of the net number of electron transfer kL through the left system-lead interface for different 
measurement times, b) The logarithm of the probability distribution highlights the behavior of the tails of the probability 
distribution. The measurement starts when the junction is at steady state and the different symbols correspond to different 
measurement times. The solid lines represent Gaussian fits. eV = 0.5 and (5 = 1. 
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FIG. 3: (a) The FT predicts that the logarithm of the ratio of the probability to measure a net number of electron transfer 
fci in the direction favored by the bias with the probability of measuring the opposite number — fc_L (which means that a net 
transfer of fcL electron occurred in the direction unfavored by the bias) is given by f3eVkL- This is given by the solid black 
diagonal line. Different symbols correspond to different measurement times (the symbols not specifically labeled correspond 
to the same measurement times as in Fig. [2}. The FT requires a sufficiently large measurement time to hold. The solid lines 
correspond to linear fits, (b) The stars represent the values of at obtained from the linear fit of the curves from (a) and the 
dashed line is the value of at predicted by Ea.p5)). 
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FIG. 4: Generating function calculated for different values of the bias and at measurement times tm such that at„ — 0.1/3eV. 
tm is calculated by solving G(tm, 0.9/3eV) = 1. The solid line is G{tm,^L) and the dotted line is G(tm, 0.9/3eV — Al). The 
difference between the two curves is a measure of the breakdown of the symmetry . 
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FIG. 5: a) Average steady state net current (/)st- The bias is at resonance with the dot levels for eV = ei = 1.697 and 
eV = £2 = 5.303. b) Average current due to the electron exiting the dots from the left side of the junction (processes = 1 on 
Fig. [!}. c) Average current due to the electron entering the dots from the left side of the junction (processes = 2 on Fig. [!}. 
d) Zero frequency power spectrum of the steady state net current [Ai + A2 using (|35|l ]. e) Mandel parameter associated to the 
process ri — 1. f) Mandel parameter associated to the process rj — 2. All quantities plotted in this figure are functions of the 
bias and are represented for five different temperature. 



